Modificamos los datos para ver las medias por artista:
spotify <- spotify %>%
select(artist, title, popularity) %>%
mutate(artist = fct_reorder(artist, popularity, .fun = 'mean'))
artist_means <- spotify %>%
group_by(artist) %>%
summarize(count = n(), popularity = mean(popularity))
Los artistas más y menos populares:
artist_means %>%
slice(1:2, 43:44)
## # A tibble: 4 × 3
## artist count popularity
## <fct> <int> <dbl>
## 1 Mia X 4 13.2
## 2 Chris Goldarg 10 16.4
## 3 Lil Skies 3 79.3
## 4 Camilo 9 81
Tenemos de 2 a 40 canciones por artista:
artist_means %>%
summarize(min(count), max(count))
## # A tibble: 1 × 2
## `min(count)` `max(count)`
## <int> <int>
## 1 2 40
La estimación no paramétrica de la densidad tiene esta pinta:
ggplot(spotify, aes(x = popularity)) +
geom_density(linewidth = 1, fill = "steelblue", alpha = .3) +
theme_bw()
Ahora vamos a ajustar el modelo full pooleado:
spotify_complete_pooled <- stan_glm(
popularity ~ 1,
data = spotify, family = gaussian,
prior_intercept = normal(50, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 84735)
Si bien en el texto dice que el prior de \(\mu\) es \(\mu
\sim N(50, 52^2)\), podemos ver que en el código dice
normal(50, 2.5, autoscale = TRUE). Esto es porque el autor
aprovecha la opción autoscale de stan_glm().
En el caso del prior gaussiano esto significa que \(\sigma\) es \(2.5
\times S_y\):
round(2.5 * sd(spotify$popularity))
## [1] 52
prior_aux por defecto si el prior de intercept es normal
es \(\sigma\). Para setear el
prior_aux también aprovecha la opción
autoscale de stan_glm(). Esto significa, para
el caso exponencial, que el prior para \(\sigma\) es \(1 /
S_y\):
round(1/sd(spotify$popularity), digits = 3)
## [1] 0.048
De hecho, si chequeamos los priors podemos ver los valore ajustados:
# Get prior specifications
prior_summary(spotify_complete_pooled)
## Priors for model 'spotify_complete_pooled'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 50, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 50, scale = 52)
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.048)
## ------
## See help('prior_summary.stanreg') for more details
Y ver el posterior summary:
complete_summary <- tidy(spotify_complete_pooled,
effects = c("fixed", "aux"),
conf.int = TRUE, conf.level = 0.80)
complete_summary
## # A tibble: 3 × 5
## term estimate std.error conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 58.4 1.10 57.0 59.8
## 2 sigma 20.7 0.776 19.7 21.7
## 3 mean_PPD 58.4 1.57 56.4 60.4
Las estimaciónes del modelo en relación a las medias por artista:
set.seed(84735)
predictions_complete <- posterior_predict(spotify_complete_pooled,
newdata = artist_means)
ppc_intervals(artist_means$popularity, yrep = predictions_complete,
prob_outer = 0.80) +
ggplot2::scale_x_continuous(labels = artist_means$artist,
breaks = 1:nrow(artist_means)) +
theme_bw() +
xaxis_text(angle = 90, hjust = 1) +
labs(x = "Artista", y = "Popularidad")
Las estimaciones promedio por artista van a ser todas iguales e iguales a las estimaciones por canción, ya que todas las estimaciones son iguales.
\[ E(\mu|y) \approx 58.39 \]
Ahora vamos a conservar todos los artistas por separado.
ggplot(spotify, aes(x = popularity, group = artist)) +
geom_density(linewidth = .5, fill = "steelblue", alpha = .3) +
theme_bw()
\[ Y_{ij}|\mu_j,\sigma \sim N(\mu, \sigma^2) \]
Ahora definimos el modelo. Lo único que cambia es que agregamos
artist como predictor (y con el -1 le sacamos
el intercept). Fijensé que ahora en lugar de
prior_intercept dice sólo prior y por eso es
que en la definici’on del libro dice que \(\mu_j \sim N(50, S^2_j)\).
spotify_no_pooled <- stan_glm(
popularity ~ artist - 1,
data = spotify, family = gaussian,
prior = normal(50, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 84735)
Acá también podemos chequear los priors:
# Get prior specifications
prior_summary(spotify_no_pooled)
## Priors for model 'spotify_no_pooled'
## ------
##
## Coefficients
## Specified prior:
## ~ normal(location = [50,50,50,...], scale = [2.5,2.5,2.5,...])
## Adjusted prior:
## ~ normal(location = [50,50,50,...], scale = [484.78,309.30,434.23,...])
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.048)
## ------
## See help('prior_summary.stanreg') for more details
Y ahora podemos ver los ajustes por artista con respecto a la \(\bar{y}_j\):
# Simulate the posterior predictive models
set.seed(84735)
predictions_no <- posterior_predict(
spotify_no_pooled, newdata = artist_means)
# Plot the posterior predictive intervals
ppc_intervals(artist_means$popularity, yrep = predictions_no,
prob_outer = 0.80) +
ggplot2::scale_x_continuous(labels = artist_means$artist,
breaks = 1:nrow(artist_means)) +
theme_bw() +
geom_point(x = 1:nrow(artist_means), y = complete_summary$estimate[1], color = "darkorange") +
xaxis_text(angle = 90, hjust = 1) +
labs(x = "Artista", y = "Popularidad")
Qué pasa si poneoms el mismo prior para \(\mu\) pero con menor dispersión. Pongamos \(0.1\) en la función que serÃa equivalente a decir que \(\mu_j \sim N(50, 2)\):
spotify_no_pooled_smallsigma <- stan_glm(
popularity ~ artist - 1,
data = spotify, family = gaussian,
prior = normal(50, .1, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
chains = 4, iter = 5000*2, seed = 84735)
# Simulate the posterior predictive models
set.seed(84735)
predictions_no <- posterior_predict(
spotify_no_pooled_smallsigma, newdata = artist_means)
# Plot the posterior predictive intervals
ppc_intervals(artist_means$popularity, yrep = predictions_no,
prob_outer = 0.80) +
ggplot2::scale_x_continuous(labels = artist_means$artist,
breaks = 1:nrow(artist_means)) +
theme_bw() +
geom_point(x = 1:nrow(artist_means), y = complete_summary$estimate[1], color = "darkorange") +
xaxis_text(angle = 90, hjust = 1) +
labs(x = "Artista", y = "Popularidad")
Efectivamente vemos que los \(\mu_j\) quedan más cerca de \(50\) y ya no son la media de las observaciones para ese artista.
Ahora tenemos básicamente tres capas de modelado:
\[ \begin{array} _Capa 1: Y_{ij}|\mu_j,\sigma_y &\sim& Modela \,\, cómo \,\, varÃa \,\, la \,\, popularidad \,\, dentro \,\, del \,\, artista_j \\ Capa 2: \mu_j|\mu,\sigma_\mu &\sim& Modela \,\, cómo \,\, varÃa \,\, la \,\, popularidad \,\, \mu_j \,\, entre \,\, artistas \\ Capa 3: \mu,\sigma_y,\sigma_\mu &\sim& Priors \,\, para \,\, los \,\, parámetros \,\, globales \,\, compartidos \end{array} \]
En la definición formal del Christensen:
\[ \begin{array} _y|\theta,\alpha &\sim& f(y|\theta,\alpha) \\ \theta|\alpha &\sim& p_0(\theta|\alpha) \\ \alpha &\sim& p_1(\alpha) \equiv p_1(\alpha|\beta_0)\\ \end{array} \]
Podemos ver la densidad de las medias por artista:
ggplot(artist_means, aes(x = popularity)) +
geom_density(linewidth = 1, fill = "steelblue", alpha = .3) +
theme_bw()
Ajustamos el modelos jerárquico. El término (1 | artist)
indica que artista es un factor de agrupamiento (o efecto aleatorio).
prior_covariance tiene que ver con la matriz de covarianza
de los efectos aleatorios, serÃa el prior de \(\mu_j\) y, no entiendo bien por qué, serÃa
lo mismo que poner un prior \(\mathcal{E}(1)\):
spotify_hierarchical <- stan_glmer(
popularity ~ (1 | artist),
data = spotify, family = gaussian,
prior_intercept = normal(50, 2.5, autoscale = TRUE),
prior_aux = exponential(1, autoscale = TRUE),
prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
chains = 4, iter = 5000*2, seed = 84735)
Y vemos los priors:
# Confirm the prior tunings
prior_summary(spotify_hierarchical)
## Priors for model 'spotify_hierarchical'
## ------
## Intercept (after predictors centered)
## Specified prior:
## ~ normal(location = 50, scale = 2.5)
## Adjusted prior:
## ~ normal(location = 50, scale = 52)
##
## Auxiliary (sigma)
## Specified prior:
## ~ exponential(rate = 1)
## Adjusted prior:
## ~ exponential(rate = 0.048)
##
## Covariance
## ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details
Podemos ver el diagnóstico sobre los 47 parámetros (44 \(\mu_j\), \(\mu\), \(\sigma_y\) y \(\sigma_\mu\)):
mcmc_trace(spotify_hierarchical) + theme_bw()
mcmc_dens_overlay(spotify_hierarchical) + theme_bw()
neff_ratio(spotify_hierarchical)
## (Intercept)
## 0.17855
## b[(Intercept) artist:Mia_X]
## 0.77280
## b[(Intercept) artist:Chris_Goldarg]
## 0.51375
## b[(Intercept) artist:Soul&Roll]
## 0.66455
## b[(Intercept) artist:Honeywagon]
## 0.79880
## b[(Intercept) artist:Röyksopp]
## 0.75600
## b[(Intercept) artist:Freestyle]
## 0.80445
## b[(Intercept) artist:DA_Image]
## 0.81100
## b[(Intercept) artist:Jean_Juan]
## 0.64785
## b[(Intercept) artist:TV_Noise]
## 0.41910
## b[(Intercept) artist:Kid_Frost]
## 0.83115
## b[(Intercept) artist:Tamar_Braxton]
## 0.84595
## b[(Intercept) artist:BUNT.]
## 0.87730
## b[(Intercept) artist:MTK]
## 0.74510
## b[(Intercept) artist:Atlas_Genius]
## 0.68040
## b[(Intercept) artist:Jazzinuf]
## 0.62940
## b[(Intercept) artist:Elisa]
## 0.79290
## b[(Intercept) artist:House_Of_Pain]
## 0.68000
## b[(Intercept) artist:Black_Stone_Cherry]
## 0.65250
## b[(Intercept) artist:C-Kan]
## 0.77330
## b[(Intercept) artist:Zeds_Dead]
## 0.68750
## b[(Intercept) artist:David_Lee_Roth]
## 0.86805
## b[(Intercept) artist:NODE]
## 0.67540
## b[(Intercept) artist:Michael_Kiwanuka]
## 0.80565
## b[(Intercept) artist:The_Wrecks]
## 0.85610
## b[(Intercept) artist:The_xx]
## 0.51225
## b[(Intercept) artist:Placebo]
## 0.67910
## b[(Intercept) artist:Mike_WiLL_Made-It]
## 0.63020
## b[(Intercept) artist:Kendrick_Lamar]
## 0.31175
## b[(Intercept) artist:X_Ambassadors]
## 0.57580
## b[(Intercept) artist:Hinder]
## 0.62940
## b[(Intercept) artist:Au/Ra]
## 0.66690
## b[(Intercept) artist:Missy_Elliott]
## 0.70295
## b[(Intercept) artist:The_Blaze]
## 0.81110
## b[(Intercept) artist:Vampire_Weekend]
## 0.69715
## b[(Intercept) artist:Alok]
## 0.37185
## b[(Intercept) artist:León_Larregui]
## 0.83365
## b[(Intercept) artist:Sufjan_Stevens]
## 0.81400
## b[(Intercept) artist:Beyoncé]
## 0.34905
## b[(Intercept) artist:Frank_Ocean]
## 0.28955
## b[(Intercept) artist:Sean_Kingston]
## 0.93845
## b[(Intercept) artist:J._Cole]
## 0.49660
## b[(Intercept) artist:Camila_Cabello]
## 0.30635
## b[(Intercept) artist:Lil_Skies]
## 0.84850
## b[(Intercept) artist:Camilo]
## 0.58210
## sigma
## 0.84050
## Sigma[artist:(Intercept),(Intercept)]
## 0.19085
rhat(spotify_hierarchical)
## (Intercept)
## 1.0008220
## b[(Intercept) artist:Mia_X]
## 1.0001391
## b[(Intercept) artist:Chris_Goldarg]
## 0.9998651
## b[(Intercept) artist:Soul&Roll]
## 1.0002406
## b[(Intercept) artist:Honeywagon]
## 1.0001168
## b[(Intercept) artist:Röyksopp]
## 1.0000496
## b[(Intercept) artist:Freestyle]
## 1.0000014
## b[(Intercept) artist:DA_Image]
## 1.0000454
## b[(Intercept) artist:Jean_Juan]
## 1.0000395
## b[(Intercept) artist:TV_Noise]
## 1.0000438
## b[(Intercept) artist:Kid_Frost]
## 1.0000608
## b[(Intercept) artist:Tamar_Braxton]
## 0.9999553
## b[(Intercept) artist:BUNT.]
## 0.9999160
## b[(Intercept) artist:MTK]
## 1.0005003
## b[(Intercept) artist:Atlas_Genius]
## 1.0001059
## b[(Intercept) artist:Jazzinuf]
## 0.9999770
## b[(Intercept) artist:Elisa]
## 0.9999609
## b[(Intercept) artist:House_Of_Pain]
## 1.0000742
## b[(Intercept) artist:Black_Stone_Cherry]
## 0.9999885
## b[(Intercept) artist:C-Kan]
## 0.9999436
## b[(Intercept) artist:Zeds_Dead]
## 0.9999352
## b[(Intercept) artist:David_Lee_Roth]
## 1.0001916
## b[(Intercept) artist:NODE]
## 1.0000100
## b[(Intercept) artist:Michael_Kiwanuka]
## 1.0000025
## b[(Intercept) artist:The_Wrecks]
## 0.9999947
## b[(Intercept) artist:The_xx]
## 1.0000460
## b[(Intercept) artist:Placebo]
## 1.0000410
## b[(Intercept) artist:Mike_WiLL_Made-It]
## 1.0000320
## b[(Intercept) artist:Kendrick_Lamar]
## 1.0004022
## b[(Intercept) artist:X_Ambassadors]
## 1.0001418
## b[(Intercept) artist:Hinder]
## 1.0002135
## b[(Intercept) artist:Au/Ra]
## 1.0000875
## b[(Intercept) artist:Missy_Elliott]
## 1.0000266
## b[(Intercept) artist:The_Blaze]
## 0.9999593
## b[(Intercept) artist:Vampire_Weekend]
## 0.9998935
## b[(Intercept) artist:Alok]
## 1.0003676
## b[(Intercept) artist:León_Larregui]
## 0.9999028
## b[(Intercept) artist:Sufjan_Stevens]
## 1.0000550
## b[(Intercept) artist:Beyoncé]
## 1.0002370
## b[(Intercept) artist:Frank_Ocean]
## 1.0006008
## b[(Intercept) artist:Sean_Kingston]
## 1.0000611
## b[(Intercept) artist:J._Cole]
## 1.0001423
## b[(Intercept) artist:Camila_Cabello]
## 1.0005154
## b[(Intercept) artist:Lil_Skies]
## 0.9999641
## b[(Intercept) artist:Camilo]
## 1.0000691
## sigma
## 1.0001962
## Sigma[artist:(Intercept),(Intercept)]
## 1.0007369
Y también la posterior:
pp_check(spotify_hierarchical) +
xlab("Popularidad") +
theme_bw()
Podemos ver también las posteriores por artista:
set.seed(84735)
predictions_hierarchical <- posterior_predict(spotify_hierarchical,
newdata = artist_means)
# Posterior predictive plots
ppc_intervals(artist_means$popularity, yrep = predictions_hierarchical,
prob_outer = 0.80) +
ggplot2::scale_x_continuous(labels = artist_means$artist,
breaks = 1:nrow(artist_means)) +
theme_bw() +
xaxis_text(angle = 90, hjust = 1) +
geom_hline(yintercept = 52.5, color = "darkorange") +
labs(x = "Artista", y = "Popularidad")
artist_summary <- tidy(spotify_hierarchical, effects = "ran_vals",
conf.int = TRUE, conf.level = 0.80)
Está interesante ver el tema de las variabilidades.
\[ \begin{array} _\frac{\sigma_y^2}{\sigma_\mu^2+\sigma_y^2} \\ \frac{\sigma_\mu^2}{\sigma_\mu^2+\sigma_y^2} \end{array} \]
artist_summary <- tidy(spotify_hierarchical, effects = "ran_vals",
conf.int = TRUE, conf.level = 0.80)
artist_chains <- spotify_hierarchical %>%
spread_draws(`(Intercept)`, b[,artist]) %>%
mutate(mu_j = `(Intercept)` + b)
artist_summary_scaled <- artist_chains %>%
select(-`(Intercept)`, -b) %>%
mean_qi(.width = 0.80) %>%
mutate(artist = fct_reorder(artist, mu_j))
ggplot(artist_summary_scaled,
aes(x = artist, y = mu_j, ymin = .lower, ymax = .upper)) +
geom_pointrange() +
theme_bw() +
xaxis_text(angle = 90, hjust = 1) +
labs(x = "Artista", y = "mu_j")
set.seed(84735)
prediction_shortcut <- posterior_predict(
spotify_hierarchical,
newdata = data.frame(artist = c("Frank Ocean", "Mohsen Beats")))
# Posterior predictive model plots
mcmc_areas(prediction_shortcut, prob = 0.8) +
ggplot2::scale_y_discrete(labels = c("Frank Ocean", "Mohsen Beats")) + theme_bw() + labs(x = "Popularidad")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.